1. Technical Field
The present invention relates to a system estimation method and program, a recording medium, and a system estimation device, and particularly to a system estimation method and program, a recording medium, and a system estimation device, in which the generation of robustness in state estimation and the optimization of a forgetting factor are simultaneously realized by using a fast H∞ filtering algorithm of a hyper H∞ filter developed on the basis of an H∞ evaluation criterion.
2. Background Art
In general, system estimation means estimating a parameter of a mathematical model (transfer function, impulse response, etc.) of an input/output relation of a system based on input/output data. Typical application examples include an echo canceller in international communication, an automatic equalizer in data communication, an echo canceller and sound field reproduction in a sound system, active noise control in a vehicle etc. and the like. For more information, see non-patent document 1: “DIGITAL SIGNAL PROCESSING HANDBOOK” 1993, The Institute of Electronics, Information and Communication Engineers, and the like.
(Basic Principle)
FIG. 8 shows an example of a structural view for system estimation (unknown system may be expressed by an IIR (Infinite Impulse Response) filter).
This system includes an unknown system 1 and an adaptive filter 2. The adaptive filter 2 includes an FIR digital filter 3 and an adaptive algorithm 4.
Hereinafter, an example of an output error method to identify the unknown system 1 will be described. Here, uk denotes an input of the unknown system 1, dk denotes an output of the system, which is a desired signal, and d^k denotes an output of the filter. (Incidentally, “^” means an estimated value and should be placed directly above a character, however, it is placed at the upper right of the character for input convenience. The same applies hereinafter.).
Since an impulse response is generally used as a parameter of an unknown system, the adaptive filter adjusts a coefficient of the FIR digital filter 3 by the adaptive algorithm so as to minimize an evaluation error ek=dk−d^k of the figure.
Besides, conventionally, a Kalman filter based on an update expression (Riccati equation) of an error covariance matrix has been widely used for the estimation of a parameter (state) of a system. The details are disclosed in non-patent document 2: S. Haykin: Adaptive filter theory, Prentice-Hall (1996) and the like.
Hereinafter, the basic principle of the Kalman filter will be described.
A minimum variance estimate x^k|k of a state xk of a linear system expressed in a state space model as indicated by the following expression:xk+1=ρ−1/2xk, yk=Hkxk+vk  (1)is obtained by using an error covariance matrix Σ^k|k−1 of the state as follows.
                                                        x              ^                                      k              |              k                                =                                                    x                ^                                            k                |                                  k                  -                  1                                                      +                                          K                k                            ⁡                              (                                                      y                    k                                    -                                                            H                      k                                        ⁢                                                                  x                        ^                                                                    k                        |                                                  k                          -                          1                                                                                                                    )                                                    ⁢                                  ⁢                                            x              ^                                                      k                +                1                            |              k                                =                                    ρ                              -                                  1                  2                                                      ⁢                                          x                ^                                            k                |                k                                                                        (        2        )                                                      K            k                    =                                                    Σ                ^                                            k                |                                  k                  -                  1                                                      ⁢                                                            H                  k                  T                                ⁡                                  (                                      ρ                    +                                                                  H                        k                                            ⁢                                                                        Σ                          ^                                                                          k                          |                                                      k                            -                            1                                                                                              ⁢                                              H                        k                        T                                                                              )                                                            -                1                                                    ⁢                                  ⁢                                            Σ              ^                                      k              |              k                                =                                                    Σ                ^                                            k                |                                  k                  -                  1                                                      -                                          K                k                            ⁢                              H                k                            ⁢                                                Σ                  ^                                                  k                  |                                      k                    -                    1                                                                                                          (        3        )                                                                    Σ              ^                                                      k                +                1                            |              k                                =                                                    Σ                ^                                            k                |                k                                      /            ρ                          ⁢                                  ⁢                  where          ,                                    (        4        )                                                                    x              ^                                      0              |                              -                1                                              =          0                ,                                  ⁢                                            Σ              ^                                      0              |                              -                1                                              =                                    ɛ              0                        ⁢            I                          ,                                  ⁢                              ɛ            0                    >          0                                    (        5        )                xk: State vector or simply a state; unknown and this is an object of estimation.    yk: Observation signal; input of a filter and known.    Hk: Observation matrix; known.    Vk: Observation noise; unknown.    ρ: Forgetting factor; generally determined by trial and error.    Kk: Filter gain; obtained from matrix Σ^k|k−1.    Σ^k|k: Corresponds to the covariance matrix of an error of x^k|k; obtained by a Riccati equation.    Σ^k+1|k: Corresponds to the covariance matrix of an error of    x^k+1|k; obtained by the Riccati equation.Σ^1|0: Corresponds to the covariance matrix in an initial state; although originally unknown, ε0I is used for convenience.
The present inventor has already proposed a system identification algorithm by a fast H∞ filter (see patent document 1). This is such that an H∞ evaluation criterion is newly determined for system identification, and a fast algorithm for the hyper H∞ filter based thereon is developed, while a fast time-varying system identification method based on this fast H∞ filtering algorithm is proposed. The fast H∞ filtering algorithm can track a time-varying system which changes rapidly with a computational complexity of O (N) per unit-time step. It matches perfectly with a fast Kalman filtering algorithm at the limit of the upper limit value. By the system identification as stated above, it is possible to realize the fast real-time identification and estimation of the time-invariant and time-varying systems.
Incidentally, with respect to methods normally known in the field of the system estimation, see, for example, non-patent documents 2 and 3.
(Applied Example to Echo Canceller)
In a long distance telephone circuit such as an international telephone, a four-wire circuit is used from the reason of signal amplification and the like. On the other hand, since a subscriber's circuit has a relatively short distance, a two-wire circuit is used.
FIG. 9 is an explanatory view concerning a communication system and an echo. A hybrid transformer as shown in the figure is introduced at a connection part between the two-wire circuit and the four-wire circuit, and impedance matching is performed. When the impedance matching is complete, a signal (sound) from a speaker B reaches only a speaker A. However, in general, it is difficult to realize the complete matching, and there occurs a phenomenon in which part of the received signal leaks to the four-wire circuit, and returns to the receiver (speaker A) after being amplified. This is an echo (echo). As a transmission distance becomes long (as a delay time becomes long), the influence of the echo becomes large, and the quality of a telephone call is remarkably deteriorated (in the pulse transmission, even in the case of short distance, the echo has a large influence on the deterioration of a telephone call).
FIG. 10 is a principle view of an echo canceller.
Then, as shown in the figure, the echo canceller (echo canceller) is introduced, an impulse response of an echo path is successively estimated by using a received signal which can be directly observed and an echo, and a pseudo-echo obtained by using it is subtracted from the actual echo to cancel the echo and to remove it.
The estimation of the impulse response of the echo path is performed so that the mean square error of a residual echo ek becomes minimum. At this time, elements to interfere with the estimation of the echo path are circuit noise and a signal (sound) from the speaker A. In general, when two speakers simultaneously start to speak (double talk), the estimation of the impulse response is suspended. Besides, since the impulse response length of the hybrid transformer is about 50 [ms], when the sampling period is made 125 [μs], the order of the impulse response of the echo path becomes actually about 400.
Non-patent document 1
“DIGITAL SIGNAL PROCESSING HANDBOOK” 1993 The Institute of Electronics, Information and Communication Engineers
Non-patent document 2
S. Haykin: Adaptive filter theory, Prentice-Hall (1996)
Non-patent document 3
B. Hassibi, A. H. Sayed, and T. Kailath: “Indefinite-Quadratic Estimation and Control”, SIAM (1996)
Patent document 1 
JP-A-2002-135171 